Modelos lineales generalizados (GLMs)

Introducción

Definición
Estructura

2.1. Función de vínculo

2.2. Familia de distribución de errores

1. Definición

Supuestos modelos lineales:

Los errores se distribuyen nromalmente
La varianza (residual) es constante
La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s), si alguna de estas variables es continua

Uno o varios supuestos no se cumplen.

¿Qué hacer?

GLMs: extensión modelos lineales especificando estructuras no normales de errores y varianzas no constantes.

Modelos lineales caso particular de GLMs con distribución de errores normal y varianza constante.

Relación varianza y media (Cayuela y de la Cruz, 2022)

Tipos de variables que por definición violan supuestos modelos lineales:

conteo de casos (p.ej. abundancia de especies)
conteo de casos como proporciones (p.ej. porcentaje plántulas muertas en un experimento)
respuesta binaria (p.ej. usado o no usado, presente o ausente)

2. Estructura

3 componentes principales:

1. El predictor lineal

2. La función de vínculo

3. La estructura de los errores

Estructura modelos lineales:

y_i = β₀ + β₁x₁ + … + β₁x_n + Ɛ_i

Ɛ_i ~N(0,σ²)

Estructura GLMs:

f(y_i) = β₀ + β₁x₁ + … + β₁x_n + Ɛ_i

Var(y) ~ f(µ)

2.1. La función de vínculo f(y)

Transformación de la variable y buscando linealizar relación y con x y/o transformando en continuas variables que no lo son.

Por tanto, cambia la naturaleza de la variable respuesta.

Determinadas funciones de vínculo deben implementarse con determinados tipos de variables respuesta

Función de vínculo (Cayuela y de la Cruz, 2022)

2.2. La familia de distribución de errores

Poisson: varianza aumenta linealmente con la media. Para respuestas tipo conteo.
Binomial negativa: extensión de Poisson. Para variables conteo con sobredispersión.
Binomial: para proporsiones y datos presencia/ausencia.
Gamma: para datos continuos con coeficiente de variación constante.

Ajuste GLM en R

Función glm()

Argumento family:

binominal(link = “logit”)
gaussian(link = “identity)
Gamma(link = “inverse”)
inverse.gaussian(link = “1/mu^2)
poisson(link = “log”)
quasi(link = “identity”, variance = “constant”)

Ejemplo con datos tipo conteo

Efecto ganadería en regeneración de araucaria. Muestreo de número de plantas de araucaria en dos tipos de ganadería: pequeños propietarios, con uso más intensivo de la tierra, y grandes empresas forestales, que mantienen algunos remanente bosque nativo. Actividad ganado medida por número de bostas.

library(ADER)
data(ara)
str(ara)

'data.frame':   26 obs. of  3 variables:
 $ seedlings: int  7 10 0 5 54 0 30 14 1 2 ...
 $ dungs    : int  5 0 330 960 115 2330 475 10 2000 130 ...
 $ property : Factor w/ 2 levels "Campesino","Empresa Forestal": 2 2 2 2 2 2 2 2 2 1 ...

Hipótesis:

H0: β = 0 (pendiente entre plántulas y número de bostas es 0)

H0: µ₁ = µ₂ (número plántulas en pequeños propietarios y empresas sin diferencias)

H0: β₁ = β₂ (relación entre plántulas y bostas es igual en pequeños propietarios y empresas)

Exploración gráfica de los datos

library(ggplot2)
ggplot(ara, aes(x = dungs, y = seedlings)) +
  geom_point(color = "steelblue", size = 2) +
  facet_wrap(~ property) +
  labs(
    title = "Relación nº plántulas y nº bostas",
    x = "Nº bostas",
    y = "Nº plántulas"
  ) +
  theme_minimal()

Construcción modelo

glm_arau1 <- glm(seedlings ~ dungs*property, ara, family = "poisson")
anova(glm_arau1)

Analysis of Deviance Table

Model: poisson, link: log

Response: seedlings

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                              25     347.56              
dungs           1   30.967        24     316.59 2.624e-08 ***
property        1   80.598        23     236.00 < 2.2e-16 ***
dungs:property  1   37.974        22     198.02 7.170e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Devianza explicada por el modelo

Medida de cantidad de variabilidad que explica el modelo (similar R² en modelos lineales)

D² = 1 - D_modelo/D_nulo

print(D2 <- 1 - (198.02/347.56))

[1] 0.4302566

summary(glm_arau1)


Call:
glm(formula = seedlings ~ dungs * property, family = "poisson", 
    data = ara)

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     2.497423   0.135408  18.444  < 2e-16 ***
dungs                          -0.010211   0.001789  -5.706 1.15e-08 ***
propertyEmpresa Forestal        0.576261   0.170467   3.380 0.000724 ***
dungs:propertyEmpresa Forestal  0.009041   0.001803   5.015 5.31e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 347.56  on 25  degrees of freedom
Residual deviance: 198.02  on 22  degrees of freedom
AIC: 282.3

Number of Fisher Scoring iterations: 5

Fórmula: log(y) = 2.497 + 0.576 Empresa + (-0.010 + 0.009 Empresa) Bostas

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_arau1)

library(aods3)
gof(glm_arau1)

D  = 198.022, df = 22, P(>D) = 2.772377e-30
X2 = 248.4275, df = 22, P(>X2) = 2.969764e-40

Modelo binomial negativo

library(MASS)
glm_arau2 <- glm.nb(seedlings ~ dungs*property, ara)
summary(glm_arau2)


Call:
glm.nb(formula = seedlings ~ dungs * property, data = ara, init.theta = 1.13715532, 
    link = log)

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     2.152593   0.362295   5.942 2.82e-09 ***
dungs                          -0.005663   0.002468  -2.294   0.0218 *  
propertyEmpresa Forestal        1.113371   0.557993   1.995   0.0460 *  
dungs:propertyEmpresa Forestal  0.003945   0.002528   1.561   0.1186    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(1.1372) family taken to be 1)

    Null deviance: 49.506  on 25  degrees of freedom
Residual deviance: 28.040  on 22  degrees of freedom
AIC: 155.22

Number of Fisher Scoring iterations: 1

              Theta:  1.137 
          Std. Err.:  0.384 

 2 x log-likelihood:  -145.218

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_arau2)

gof(glm_arau2)

D  = 28.0403, df = 22, P(>D) = 0.1743489
X2 = 24.0659, df = 22, P(>X2) = 0.343783

Predicciones modelo

library(car)
Anova(glm_arau2, type = "III")

Analysis of Deviance Table (Type III tests)

Response: seedlings
               LR Chisq Df Pr(>Chisq)   
dungs            7.9871  1   0.004711 **
property         4.3145  1   0.037789 * 
dungs:property   3.5106  1   0.060977 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(visreg)
visreg(glm_arau2, "dungs", by ="property", overlay = T, scale = "response", ylab = "Nº plántulas", xlab = "Nº bostas")

Ejemplo con distribución de errores binomial

2 variantes de estos modelos:

Variable con valores binarios (1/0): modelos de regresión logística
variables con valores porcentuales (nº éxitos / nº de ensayos)

Función de vínculo: logit log(p/(1-p))

logit(y) = ln(p/(1-p)) = β₀ + β₁x₁ + … + β₁x_n

Curva sigmoidea entre 0-1

Ejemplo: Estrategia reproductora especie de alga invasora

Estrategia reproductora de un alga invasora con dos especies nativas.

¿Probabilidad de encontrar estructuras reproductoras es diferente entre la especie invasora y las nativas?

Hipótesis: la especie invasora tendrá mayor probabilidad de reproducrise que las nativas

library(ADER)
data(algas)
str(algas)

'data.frame':   359 obs. of  3 variables:
 $ Sp    : Factor w/ 3 levels "Fra","Tom","Ver": 3 3 3 3 3 3 3 3 3 3 ...
 $ Long  : num  1 2 4 4 4 4.5 4.5 4.5 5 5 ...
 $ Estado: int  0 0 0 0 0 0 0 0 0 0 ...

Variable respuesta (Estado):

 - 0: sin estructuras reproductoras
 - 1: con estructuras reproductoras

Variables predictoras:

 - talla del individuo (cm)
 - especie (Sp)

Exploración de los datos

ggplot(algas, aes(x = Long, y = Estado, color = Sp)) +
  geom_point(alpha = 0.4, position = position_jitter(height = 0.05)) +
facet_wrap(~Sp) +
  labs(
    x = "Talla individuo (cm)",
    y = "Estructura reproductora",
  ) +
  theme_minimal()

Ajuste del modelo

glm_algas1 <- glm(Estado ~ Long*Sp, algas, family = "binomial")
anova(glm_algas1)

Analysis of Deviance Table

Model: binomial, link: logit

Response: Estado

Terms added sequentially (first to last)

        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                      358     476.38              
Long     1  103.010       357     373.38 < 2.2e-16 ***
Sp       2   31.771       355     341.60 1.262e-07 ***
Long:Sp  2    0.459       353     341.15    0.7949    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

glm_algas2 <- glm(Estado ~ Long+Sp, algas, family = "binomial")
anova(glm_algas2)

Analysis of Deviance Table

Model: binomial, link: logit

Response: Estado

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   358     476.38              
Long  1  103.010       357     373.38 < 2.2e-16 ***
Sp    2   31.771       355     341.60 1.262e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(D2 <- 1 - (341.6/476.38))

[1] 0.2829254

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_algas2)

Estructura del modelo

summary(glm_algas2)


Call:
glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.80703    0.42701  -6.574 4.91e-11 ***
Long         0.26505    0.03174   8.351  < 2e-16 ***
SpTom       -0.57476    0.34528  -1.665    0.096 .  
SpVer       -1.84877    0.36933  -5.006 5.57e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 476.38  on 358  degrees of freedom
Residual deviance: 341.60  on 355  degrees of freedom
AIC: 349.6

Number of Fisher Scoring iterations: 5

Gráfico de predicciones

newdat <- expand.grid(
  Long = seq(min(algas$Long), max(algas$Long), length = 100),
  Sp = levels(algas$Sp)
)

# Predicciones de probabilidad
pred <- predict(glm_algas2, newdat, type = "link", se.fit = TRUE)

# Convertir de escala logit a probabilidad
newdat$fit <- plogis(pred$fit)
newdat$lwr <- plogis(pred$fit - 1.96 * pred$se.fit)
newdat$upr <- plogis(pred$fit + 1.96 * pred$se.fit)

ggplot(newdat, aes(x = Long, y = fit, color = Sp, fill = Sp)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.15, color = NA) +
  labs(
    x = "Talla individuo (cm)",
    y = "Predicción probabilidad estructuras reproductoras",
  ) +
  theme_minimal(base_size = 13)

Comparaciones entre especies

library(multcomp)
summary(glht(glm_algas2, linfct = mcp(Sp = "Tukey")))


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)

Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)    
Tom - Fra == 0  -0.5748     0.3453  -1.665 0.217801    
Ver - Fra == 0  -1.8488     0.3693  -5.006  < 1e-04 ***
Ver - Tom == 0  -1.2740     0.3126  -4.076 0.000138 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)